Devonshire-Landau free energy of BaTiOs from first principles 
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We have studied the Devonshire-Landau potential underlying the phase transition sequence of 
BaTiOa using the first-principles effective Hamiltonian of Zhong, Vanderbilt, and Rabe [Phys. Rev. 
Lett. 73, 1861 (1994)], which has been very successful in reproducing the phase transitions and the 
dielectric and piezoelectric properties of this compound. The configuration space (determined by 
the polarization P as order parameter) was explored with the help of auxiliary electric fields. We 
show that the typically assumed form of the potential, a sixth-order expansion in P around the 
paraelectric cubic phase, properly accounts for the behavior of the system, but we find a non-trivial 
temperature dependence for all the coefficients in the expansion, including the quadratic one, which 
is shown to behave non-linearly. Our results also prove that the sixth-order terms in the free-energy 
expansion (needed to account for the first-order character of the transitions and the occurrence of 
an orthorhombic phase) emerge from an interaction model that only includes terms up to fourth 
order. 



I. INTRODUCTION 

From a phcnomenological point of view, the behav- 
ior of a system in the vicinity of a phase transitionmcan 
be described in the framework of Landau theory.^ In 
this scheme, one begins by identifying the so-called or- 
der parameter, a (in general multidimensional) variable 
Q that characterizes the symmetry change in the tran- 
sition, and then constructs the Landau free-energy func- 
tion T(Q;T), with the property that the equilibrium 
value of Q as a function of temperature is that which 
minimizes T. Formally, the Landau free energy is an 
incomplete thermodynamic potential, which in principle 
could be calculated as: 



F(Q;T) 



jTln]Te-^/ fcBT 
j/Q 



(1) 



in analogy with the procedure to obtain the standard free 
energy in terms of the partition function. Here the sum 
is only over those states with the given value of Q. In 
a Landau phenomenological treatment, one just assumes 
that T(Q; T) exists, and that it can be represented as a 
simple power series in Q in the vicinity of the transition: 



F(Q;T) = F (T) + A(T)Q 2 



pO) 



•s£4 4) (Q) 



(2) 



where fg are invariants of order j constructed from Q. 
The quadratic coefficient A is assumed to vary linearly 
with temperature in the form A(T) = a(T — To), in such 
a way that To is the transition temperature (for a second- 
order transition) or the lower metastability limit of the 
upper phase (for a first-order one). These assumptions 
are known to break down in the critical region close to the 
transition point of continuous phase transformations, but 
they are valid when describing the approximate behav- 
ior of the system in wider temperature intervals around 
the transition temperature. Classical Landau theory ig- 
nores any temperature variation of the higher-order co- 



efficients in the expansion. Nevertheless, it has been ap- 
plied quite successfully to the analysis of many materials 
in wide temperature ranges.B Nowadays Landau theory is 
the most common approach to study the phenomenology 
of any symmetry-breaking structural phase transition. 

One of the early successes of this phenomenological ap- 
proach to phase transitions was the description by Devon- 
shiretl of the sequence of transitions in barium titanate 
(BaTiOa), a sequence that spans a temperature inter- 
val of more than 200 degrees. This material exhibits at 
high temperatures a paraelectric cubic perovskite struc- 
ture and, as temperature decreases, it undergoes three 
successive first-order transitions to ferroelectric phases 
with tetragonal, orthorhombic, and rhombohedral sym- 
metries. In these phases the polarization P points along 
one of the (1,0,0), (1,1,0), and (1,1,1) cubic direc- 
tions respectively^ The polarization P can be identi- 
fied as the (three-dimensional) order parameter, and the 
Devonshire-Landau potential (per unit volume) written 
as: 



F = 



T + \aP 2 + uP A 
+h 1 P 6 + h 2 (P^ + 
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where To is the free energy of the reference cubic phase 
and P 2 stands for P 2 + P 2 + P 2 (the notation is taken 
from the review paper by Cowleju). This is the com- 
plete expansion of T up to sixth order, and its relatively 
simple form is due to the high symmetry of the system. 
Sixth-order terms in P in the expansion are needed to 
account for the first-order character of the transitions 
(to model non-equivalent coexisting free-energy minima). 
Moreover, it can be seen that one needs anisotropic sixth- 
order terms, in T to account for the orthorhombic phase 
of BaTi03.Q Devonshire showed that by assuming a lin- 
ear temperature dependence for a and suitable constant 
values for the rest of the coefficients it is possible to qual- 
itatively reproduce the transition sequence as well as the 
dielectric properties of the system. 



After Devonshire's work, the use of an expansion like 
that of Eq. y has become the standard way to model the 
thermodynamic properties of ferroelectric perovskites, 
even though the correctness of the assumption of con- 
stant high-order coefficients was questioned early on. 
Drougard et al. and Huibregtse et alB studied BaTi03 
at temperatures just above the cubic to tetragonal tran- 
sition and in the orthorhombic phase respectively; they 
assumed the existence of a Devonshire-Landau potential 
but showed that there should be a significant tempera- 
ture dependence of the fourth-order termsE It has since 
been shown that there is a relatively large latitude to pro- 
duce Devonshire-like models that lead to qualitatively 
sensible predictions for the transition sequence, diver- 
gence of the dielectric constants, etc., even though these 
models may not be realistic when confronted with other 
experimental evidence£3 On the other hand, it could be 
questioned whether such a relatively simple free-energy 
expansion as that in Eq. |3| exists at all. For instance, 
Nakamura and Kinase have worked on a different ap- 
proach to the problem in which a non-polynomical from 
of the free energy is used; they also discuss the connection 
of their work to the Devonshire theoryJl^l 

In the past decade, first-principles methods based 
on density-functional theory have demonstrated to be 
accurate enough to reproduce and predict structural 
phase transitions in perovskite oxidesJia In particular, 
the phas^ptransition sequence and other properties of 
BaTi0 3 ,yy PbTi0 3 ,N and KNboJla have been stud- 
ied using the so called effective Hamiltonian approach. 
An ab-initio effective Hamiltonian H e s is a mechanical 
model that includes the relevant microscopic degrees of 
freedom of the system and is constructed on the basis 
of first-principles calculations. This model can then be 
analyzed by statistical methods (typically, Monte Carlo 
or molecular dynamics) to explore the finite-temperature 
behavior of the system. Monte-Carlo simulations with an 
ab-initio effective Hamiltonian for BaTi03 have been ex- 
tremely successful, replicatiag approximately the experi- 
mental transition sequencecJ and succeeding in reproduc- 
ing the main features of the [dielectric and piezoelectric 
properties of the real systemJlJ 

A natural question to ask then is to what extent the 
thermal behavior resulting from these ab-initio effective 
Hamiltonians is compatible with the Devonshire-Landau 
framework. In this paper we show that a Devonshire- 
Landau expansion of the form of Eq. g does indeed 
emerge from detailed Monte Carlo simulations with the 
ab-initio effective Hamiltonian for BaTi03 of Ref. O, 
with the important qualification that the coefficients 
in the expansion have a non-trivial temperature depen- 
dence. In particular, we show that the sixth-order terms 
in F(P] T) are non-zero only in the temperature range in 
which the transitions occur. Moreover, it becomes clear 
that these sixth-order coefficients appear as products of 
the statistical fluctuations of a Hamiltonian that only in- 
cludes terms up to fourth order in the polar degrees of 
freedom. 



II. METHOD AND TECHNICAL DETAILS 

Zhong, Vanderbilt, and Rabsi3 constructed their ef- 
fective Hamiltonian for BaTi03 retaining as relevant de- 
grees of freedom of the system a local polar distortion u^ 
in each cell, the homogeneous strain 77, and an inhomoge- 
neous strain represented by a second set of local vectors 
v;. The H c s has the form: 



H, 



err 



+ £i i „, /J r(a,/3)u? a «? /9 

+ Y,i- a ,p;l B ( i: > a > # l)Ui a Uipr)l 
+ T l i,hB(l,h)r 1 i7 lh , 



(4) 



where, for clarity, we have not written the terms asso- 
ciated to the inhomogeneous strain. Here i and j range 
over the cells in the system, a and j3 are cartesian indexes, 
and I and h refer to tensor components in Voigt nota- 
tion. For the local polar modes, the H e g contains har- 
monic couplings J(i,j;a,/3) (on-site, i = j, and between 
modes in different cells up to the third nearest neigh- 
bors) that xeproduce the instabilities of the cubic phase 
of BaTi03,EH and fourth-order on-site terms T(a, 0) that 
define the low-symmetry minima and stabilize the sys- 
tem. The effect of strain is included through the stan- 
dard elastic energy and through the coupling coefficients 
B(i; a, /3] I), which account for the spontaneous strain in 
the low-symmetry phases. This model is probably the 
simplest one that captures the essential physics of the 
system. As mentioned in the introduction, H B g is fourth- 
order in the polar variables. 

A direct link between this iJ e ff and the phenomeno- 
logical potential of Eq. g can be established in the very 
low-temperature limit T — ► 0. Eq. [I] shows that F(P; 0) 
is just the energy of the lowest-lying state with polariza- 
tion P. For BaTiC>3 this state exhibits homogeneously 
polarized cells with the global strain adjusted so as to 
minimize the energy (the inhomogeneous strain is zero). 
It can be shownliS that for a given P, the homogenous 
strain is proportional to the square of the polarization, so 
F(P; 0) is fourth-order in the polarization. This moans 
that h± = h% = /13 = in the low-temperature limit. Us- 
ing the H c ?i parameters, we obtain the rest of coefficients 
in F at K: a = -4/3, u = 0.094, and v = 0.051. Here, 
units have been chosen so that the first-octant rhombo- 
hedral ground state minimum is located at P = (1,1,1) 
and its energy per unit volume with respect to the cubic 
phase is —1. (This election fixes the value of a and the 
relation 9u + 3v = 1; i.e., F(P; 0) is determined only by 
one coefficient.) We have plotted this K free-energy 
map in Fig. Pi), which shows the cubic phase as a local 
free-energy maximum, tetragonal and orthorhombic sad- 
dle points (both stable in the radial direction), and the 
rhombohedral global minima. 
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FIG. 1. Free-energy maps of BaTiOa at several tempera- 
tures. Panel a): 300 K, b): 290 K, c): 250 K, d): 220 K, 
e): 210 K, f): 190 K, g): 100 K, and h): OK. For each 
temperature we show the (010) and (110) planes of the or- 
der-parameter configuration space, with contour lines de- 
picted only in the low free-energy regions. Symbols charac- 
terize the relative stability of the critical points with respect 
to displacements within a given plane: circles for free-energy 
minima (a bigger circle corresponds to the stable-equilibrium 
state), squares for free-energy maxima, and triangles for sad- 
dle points ("up" (resp. "down") triangles for points that are 
maxima (resp. minima) along the radial direction) . Note that 
critical points with orthorhombic symmetry may need to be 
represented by two different symbols in the (010) and (110) 
planes (see panel b)). 

Now we tackle the calculation of F at finite tempera- 
tures. A straightforward approach to obtain information 
about— the Landau potential from the results of Monte 
Carlollj simulations was presented by Radescu et aln in 
the context of a three dimensional $ 4 model. From Eq. |l| 
it is easy to show that there exists a relation between the 
probability distribution of the order parameter V(Q:T) 
and the corresponding Landau potential F(Q;T): 



F(Q;T) - F eq {T) = -k B Tln[V(Q;T)}, 



(5) 



where F eq (T) is the equilibrium free energy of the system 
at the temperature considered. This procedure, which 
can be quite efficiently applied to relatively simple cases 
such as the $ 4 model (one-dimensional order parameter, 
a single second-order phase transition), is not well suited 
for our problem. We have to study a three-dimensional 
free-energy map with coexisting non-equivalent minima. 
A typical Monte Carlo simulation gives information only 
about one free-energy minimum (usually, the one corre- 
sponding to the stable-equilibrium state of the system) , 
with a poor sampling of the order-parameter probability- 
distribution function in other regions. On the other hand, 
more sophisticated sampling strategies would be very de- 



manding from the computational point of view. 

This problem can be overcome, and F(P\T) computed 
in an efficient manner, by modifying the effective Hamil- 
tonian sa as to include the effect of an external electric 
field E:N 



H 'ctf 



H, 



E-Z* 



Z^ 



(6) 



where Z* is the Born effective charge associated to the 
local polar distortions. As the quantity multiplying the 
electric field is just the total dipole moment of the sys- 
tem, the new Landau potential is simply F'(P;T, E) = 
F(P; T) — EP. An applied electric field changes the 
location (and possibly the symmetry) of the stable- 
equilibrium state. The new location can be determined 
by a MC simulation of the modified effective Hamiltonian 
(P e9 = (P)) and also computed directly from F 1 . The 
idea can be mathematically expressed as: 



dF'(P; T, E) 
dP 



(7) 



Pe,(T,E) 



The left hand side of this equation depends linearly on 
the coefficients of F and on the electric field. Thus, for a 
given temperature we can consider several electric fields, 
perform MC runs to obtain the equilibrium values P eq ,E£l 
and then find the best solution of an overdetermined set 
of linear equations of the form of Eq. to get the coeffi- 
cients of the expansion of F. 

For each temperature, we first performed a zero-field 
calculation to obtain a snapshot of the stable-equilibrium 
configuration. We then used this configuration as the 
starting point for the runs with electric field applied. In 
order to illustrate the kind of information we were pur- 
suing, consider T such that the stable-equilibrium state 
of the system is orthorhombic with V eq — P eq (l, 1,0). In 
this case, we would first use fields of the form (E, E, 0) 
to obtain information about the response of the sys- 
tem within the orthorhombic phase. Fields of the forms 
(E, 0, 0), (0, 0, E), and (E, E, E) would probe the relative 
stability of minima with different symmetries. Finally, 
general fields leading to P xm ^ P y ^ eq ^ P zm ^ P x , eq , 
would explore in detail the intermediate regions of the 
configuration space. In our Monte Carlo simulations we 
used a 14 x 14 x 14 supercell (which corresponds to 13720 
atoms) with periodic boundary conditions; we typically 
performed 15000 MC sweeps for the thermalization of 
the system and 35000 more to calculate the averages of 
the polarization (these are very well-converged calcula- 
tion conditions, as can be checked in Refs. 13|-^5|). For 
each temperature we considered around one hundred dif- 
ferent electric fields and constructed a largely overdeter- 
mined system of equations for the coefficients in F, which 
could then be reliably fitted. 



III. RESULTS AND DISCUSSION 

The data from our Monte Carlo simulations can in- 
deed be represented by a Landau free energy in the form 
of Eq. in a wide temperature interval that includes all 
the transitions. (The transition temperatures predicted 
by the effective Hamiltonian are T c \ = 297 K (cubic to 
tetragonal), T C 2 — 230 K (tetragonal to orthorhombic) 
and T C 3 =j— 200 K (orthorhombic to rhombohedral) re- 
spectively.Eil) Our results show that the coefficients in the 
expansion of F have a significant and non-trivial temper- 
ature dependence. 

Before discussing in detail the behavior of the coeffi- 
cients, it is helpful to present the shape of the free energy 
they determine at a few temperatures in the range of the 
phase transitions. Fig. Ilia,) (300 K) shows how tetragonal 
local minima appear just above the first transition, while 
the absolute minimum of F is is still located at P=0 in 
the cubic well. Panels b) and c) correspond to 290 K and 
250 K respectively, i.e., temperatures in the range of sta- 
bility of the tetragonal phase. We see that the tetragonal 
wells are indeed the global free-energy minima, and that 
the orthorhombic and rhombohedral wells nucleate in the 
form of saddle points that are unstable with respect to 
a tetragonal distortion (in the rhombohedral case, the 
saddle points are unstable with respect to an orthorhom- 
bic distortion also). Panels d) and e) refer to 220 K and 
210 K respectively, i.e., temperatures in the orthorhom- 
bic range. The tetragonal wells become metastable with 
respect to the orthorhombic ones, which are now the free- 
energy global minima. At the same time, rhombohedral 
local minima appear. Finally, panels f ) and g) correspond 
to 190 K and 100 K respectively, both in the rhombohe- 
dral temperature range. Here we see that the rhombo- 
hedral minima have finally become the deepest ones. It 
is clear that the sequence converges to the free-energy 
map given by the effective Hamiltonian itself (panel h)). 
For all the transitions, the coexistence of different free- 
energy minima in the figures evidences their first-order 
character. 

Our method to explore the F landscape has a limita- 
tion when probing the vicinity of the cubic phase. After 
the point P=0 becomes a local free-energy maximum a 
few degrees below T c \, it is no longer possible to sample 
the region around it using auxiliary electric fields: the cu- 
bic structure is already unstable and thus no longer useful 
as a starting point for the simulations, and attempts to 
steer a system in a tetragonal, orthorhombic, or rhom- 
bohedral well towards P=0 only succeed in landing it 
in the corresponding symmetry-inverted domain. As the 
behavior of F around P=0 is basically determined by 
the quadratic coefficient a of the Devonshire-Landau ex- 
pansion, the value assigned to this coefficient by a fit of 
the simulation data should be considered suspect. (This 
problem pertains only to the region near P=0. Informa- 
tion about the relative stability of tetragonal, orthorhom- 
bic, and rhombohedral phases is available at any temper- 



ature, since it is always possible to apply fields that lead 
the system to each of these symmetries.) 
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FIG. 2. Temperature behavior of the quadratic coefficient 
(a) of the Devonshire-Landau potential: The dashed line con- 
nects the values obtained from a direct fit to MC data, the 
solid line shows the interpolated o(T), and the solid circles 
are the points used to fit the interpolating polynomial (see 
text). The transition temperatures are marked with vertical 
lines. 

As shown in Fig. 0, at temperatures around and above 
T c \ the behavior of a is fairly reasonable. It is posi- 
tive at high temperatures, becoming negative a few de- 
grees below T c i, as corresponds to a first-order transi- 
tion. For very low temperatures a also behaves well, 
tending smoothly to its K value. This is because in 
this region the underlying Landau potential is so sim- 
ple (the thermal fluctuations are relatively small and the 
sixth-order terms almost negligible) that we can use the 
information obtained around the low-symmetry minima 
to reconstruct the complete free-energy map. However, 
in the intermediate region (from 150 K to 250 K), a 
turns positive again, which would mean that the cubic 
phase becomes metastable in a wide temperature inter- 
val. This metastability is an unphysical result (explicitly 
ruled out within our model by zero-field MC simulations 
starting from a cubic phase, in which we found no trace of 
it). We thus reconsidered the fitting of our data, impos- 
ing a more physical temperature evolution of a: we as- 
sumed that in the intermediate-temperature range a(T) 
is given by a simple interpolation (Fig. ||) between the 
high-temperature region, where we can reliably sample 
it, and the low-temperature limit determined directly by 
the convergence to the effective Hamiltonian (a third- 
order polynomial suffices for our purposes). With this 
from for a(T) fixed, our MC data are still well fitted to 
the Landau potential of Eq. || For instance, at 190 K the 
relative error obtained when a(T) is freely fitted is 2.2%, 
and when a(T) is fixed by the interpolation the error is 
still very small: 3.6%. This clearly shows that, in the 
intermediate-temperature region, our Monte Carlo data 
contain little information about the value of a£3 Higher- 



order terms were tried but found to play no role in the 
fit, so they can be excluded from the potential. Our final 
results are shown in Fig. |[ 
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FIG. 3. Free-energy coefficients fitted to MC data assum- 
ing a(T) is given by the solid line in Fig. 0. The transition 
temperatures are marked with vertical lines. 

At very low temperatures a, u, and v tend to their 
K values and the sixth-order terms go to zero. On the 
other hand, at temperatures well above the first transi- 
tion (T > 350 K) we find the MC data are properly fit- 
ted to a fourth-order Landau potential, so the sixth-order 
terms can again be excluded from the model. The sixth- 
order terms then occur only in the temperature range 
in which the transitions take place (as plainly shown in 
Fig. 0) . It is also remarkable that the fourth-order terms 
exhibit a very strong temperature dependence. Partic- 
ularly, the anisotropic term v changes sign at around 
90 K and takes large negative values all through the 
intermediate-temperature range. So, apart from the in- 
fluence it has, together with hi and /13, in determin- 
ing the transition sequence, we find that a negative v 
is responsible for the first-order character of the cubic to 
tetragonal transition. 

In Fig. [| we have plotted our coefficients near the cubic 
to tetragonal transition temperature together with the 
available experimental results.™ (We have changed the 
temperature scale to make the experimental T c \ coincide 
with the theoretical one.) The agreement is only quali- 
tative but as good as could be expected given the sim- 
plifying assumptions involved in the experimental work 
of Ref. H (the coefficients are restricted to be constant or 
to depend linearly on temperature) and the fact that the 
differences among the two experiments are comparable 
to those between experiment and theory. 

Our work strongly suggests that the phase transitions 
of BaTiOa can be described in terms of a single Lan- 
dau potential F with the form of Eq. 0. Despite what 
is assumed in most of the previous work on this prob- 
lem, the quadratic parameter a is found to exhibit a 
strongly non-linear temperature dependence. This is 
clear because we can reliably calculate a at high (over 



T c i) and very low temperatures, and the two regions can- 
not be joined linearly. We show that all the high-order 
terms in F present a very significant evolution with tem- 
perature. This conclusion is essentially opposed to the 
temperature-independent behavior that is still assumed 
by some authors (see, for example, Ref. [Kj). It is also 
very remarkable that the two features of BaTi03 that 
require the inclusion of sixth-order terms in the Lan- 
dau potential, i.e., the first-order character of the transi- 
tions and the occurrence of an orthorhombic phase, have 
been reproduced using a fourth-order effective Hamilto- 



nian.cij This piece of information should be taken into 
account when constructing mechanical models to study 
the finite-temperature behavior of similar compounds. 
For instance, the authors of Ref. [n] included sixth-order 
terms in the underlying interaction model for BaTiC>3 
with the explicit aim to account for the first-order char- 
acter of the transitions; our results indicate that those 
terms should be unnecessary. 



a, which is shown to behave non-linearly. Our work also 
shows that the sixth-order terms in polarization needed 
to explain basic features of BaTiC>3 in a Devonshire- 
Landau approach (the first-order character of the transi- 
tions and the occurrence of an orthorhombic phase) are 
properly accounted for by an interaction model that only 
includes terms up to fourth order. 
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FIG. 4. Comparison of our free-energy coefficients (lines 
with the experimental values of Refs. M (open symbols) and 
(solid symbols): solid line and circles for a, dotted line and 
squares for u + v, dashed line and "up" triangles for hi + hz, 
and dash-dotted line and "down" triangles for h^. The verti- 
cal line marks the cubic-tetragonal transition temperature. 



IV. SUMMARY 

We have studied the Devonshire-Landau potential un- 
derlying the phase transition sequence of BaTiC>3 using 
the first-principles effective Hamiltonian parametrized by 
Zhong, Vanderbilt, and Rabe. The order-parameter con- 
figuration space was explored with the help of auxiliary 
electric fields that change the location and relative sta- 
bility of the free-energy minima. Our results show that 
the typically assumed form of the potential, an expansion 
up to sixth order in the polarization from the paraelec- 
tric cubic phase, properly accounts for the behavior of 
the system. But, despite what is usually presumed, we 
find a non-trivial temperature dependence for all the co- 
efficients in the expansion, including the quadratic term 



* Email address: wdbingoj@lg.ehu.es 

L. D. Landau and E. M. Lifshitz, Statistical Physics (Perg- 
amon Press Ltd., London, 1958). 

2 S. Radescu, I. Etxebarria, and J. M. Perez-Mato, J. Phys.: 
Cond. Matt. 7, 585 (1995). 

3 A. F. Devonshire, Phyl. Mag. 40, 1040 (1949); ibid. 42, 
1065 (1951); Advances in Physics 3, 85 (1954). 

4 See, for example: F. Jona and G. Shirane, Ferroelectric 
Crystals (Pergamon Press, Oxford, England 1962; reedi- 
tion by Dover Publications, Inc., New York 1993) or M. 
E. Lines and A. M. Glass, Principles and Applications of 
Ferroelectric Materials (Clarendon Press., Oxford 1977). 

5 R. A. Cowley, Advances in Physics 29, 1 (1980). 

6 For a discussion of the conditions that a Devonshire- 
Landau potential must fulfill to account for phase s of a 



given symmet ry, see D. Vanderbilt and M. H. Cohen, cond 



mat/0009337j . 

In the second article in Ref. pi Devonshire extended the 
theory to consider the strain and gave a proper account for 
the elastic and piezoelectric properties of the system. 

8 M. E. Drougard, R. Landauer, and D. R. Young, Phys. 
Rev. 98, 1010 (1955); E. J. Huibregtse and D. R. Young, 
ibib. 98, 1562 (1956). 

9 A temperature dependence of the sixth-order terms in F 
was reported by J. A. Gonzalo and J. M. Rivera in Ferro- 
electrics 2, 31 (1971). 

10 See the discussion in Ref. ]a and that by K. Fujita and Y. 
Ishibashi, Jpn. J. Appl. Phys. 36, 254 (1997). 

11 K. Nakamura and W. Kinase, J. Phys. Soc. Jap. 61, 4596 
(1992); ibid. 61, 2114 (1992). 

12 See the short review by D. Vanderbilt in Current Opinions 
in Solid State and Materials Science 2, 701 (1997) and ref- 
erences therein. 

13 W. Zhong, D. Vanderbilt, and K. M. Rabe, Phys. Rev. Lett. 
73, 1861 (1994); Phys. Rev. B52, 6301 (1995). 

14 A. Garcia and D. Vanderbilt, Appl. Phys. Lett. 72, 2981 
(1998); Proceeding of the First Principles Calculations for 



Ferroelectrics: Fifth Williamsburg Workshop, edited by R. 
E. Cohen (AIP, Woodbury, New York, 1998), p. 53. 

15 U. V. Waghmare and K. M. Rabe, Phys. Rev. B52, 13236 
(1995). 

16 H. Krakauer, R. Yu, C. Z. Wang, and C. Lasota, in Proceed- 
ings of the 1997 Williamsburg Workshop on Ferroelectrics, 
edited by H. Chen, Ferroelectrics 206, 133 (1998). 

17 The dispersion curves associated to the harmonic effective 
system are essentially those of the unstable band calcu- 
lated by Ghosez et ed. [Ph. Ghosez, E. Cockayne, U. V. 
Waghmare, and K. M. Rabe, Phys. Rev. B60, 836 (1999)]. 

18 R. D. King-Smith and D. Vanderbilt, Phys. Rev. B49 5828 
(1994). 

Although our emphasis is on Monte Carlo simulations, the 
discussion applies also to molecular dynamics methods. 

20 From a MC simulation we typically obtain (P) correspond- 
ing to the stable-equilibrium state of the system. Eq. |7J 
can be used also if the system gets stuck in a metastable- 
equilibrium state, but not if the simulation alternates from 
one free-energy minimum to another, giving a spurious (P). 

21 The experimental transition temperatures for BaTiOa are 
Td = 403 K, T c2 = 278 K, and T c3 = 183 K. 

22 We tried to compute a in the problematic temperature in- 
terval by modifying the effective Hamiltonian to keep the 
point P = as a local minimum (what could be under- 
stood as a kind of umbrella sampling, see M. P. Allen and 
D. J. Tildesley, Computer Simulation of Liquids (Oxford, 
New York, 1990)). However, this modification results in 
inhomogeneous equilibrium states, which lie outside the 
Devonshire-Landau framework. 

23 For Pb(Zr a: Tii_ a; )03, a similar fourth-order effective 
Hamiltonian reproduces a monoclinic phase, and thus 
would formally lead to eigth-order terms in the Devonshire- 
Landau potential. See Ref. o and references therein. 



